clc
clear

L = 5;
dx = 0.1;

[x y] = meshgrid(0:dx:L);
n = size(x, 1);

%A = zeros(n,n);
A = spalloc(n^2, n^2, 5*n^2); %根据AI建议，使用稀疏矩阵优化内存
b = zeros(n^2,1);

boundary_top = @(x,y) sin(pi/L*x); %上边界

for i = 1:n
  for j = 1:n
    _k = i+(j-1)*n;

    if i == 1 || i == n || j == 1 || j == n %如果是边界
      b(_k) = 0;
      A(_k,_k) = 1;
      if i == n %上边界
        b(_k) = boundary_top(x(i,j),y(i,j));
      end
    else
      A(_k,_k) = -4;
      A(_k,_k+1) = 1;
      A(_k,_k-1) = 1;
      A(_k,_k+n) = 1;
      A(_k,_k-n) = 1;
    end
  end
end

%_u = zeros(n^2,1);
_u = A\b; % A是稀疏阵
u = reshape(_u,n,n); %将u变回n*n的场

% 计算梯度场
grad_x = zeros(n);
grad_y = zeros(n);
grad_x(2:n-1,2:n-1) = -(u(2:n-1,3:n) - u(2:n-1,1:n-2))/(2*dx);
grad_y(2:n-1,2:n-1) = -(u(3:n,2:n-1) - u(1:n-2, 2:n-1))/(2*dx);


figure();
hold on
axis equal
h = surf(x,y,u);
set(h,'edgecolor','none');
quiver3(x,y,5*ones(n),grad_x,grad_y,zeros(n));
xlabel('x');
ylabel('y');

